In [165]:
import numpy 
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

%matplotlib inline
from pylab import *

Liddriven Cavity :

Before exploring this example note that we changed the basic code and removed loops . This really helps us to save the time . If you run this code with considering loop instead of indexes it will take approximately 15 hours to run .

In this simulation we have a 2D fluid flow that is driven by a lid at the top which moves at a speed of uo . the cavity has four sides which three of them are solid and no-slip boundary condition is considered for walls and the top one which has moving velocity . If you remember in the file (illustration) we mention how to deal with moving or solid boundaries . This problem is twodimensional .

We need to consider both advection and diffusion roles in Navier-stokes to ideally simulate this problem . In this example in contrast to previous examples the fluid has the ability to move .

Question : How can we measure the strength of motion or strength of the fluid induced by moving top lid ? We have the answer thanks to the fluid mechanic course . The fluid could be Water or Gas with different kinematic velocities or densities .

We at first need to non-dimensionalize the navier stokes equation . If we do the Reynolds number will appear which is defined as the inertia force over viscosity force . Higher reynolds number implies higher inertia or lower viscosity resistance and vice verca .

What changes should be considered to shift our code from conduction to fully Navier-Stokes equation ?

1- define parameter reynolds which you want to consider , It means that you definitely understand that you can not directly the parameters which you directly use in FDM or in original format of your equation . For example if you set the reynolds number with adjusting the top velocity or kinematic viscosity in original equation here you are not allowed to consider any arbitary value .

The most important point is that you must make equal value of Reynold number . You need to have reynold number for example 1000 then you choose the charectristic velocity as 0.1 ( large value breaks low mach number ) and also you need to have relaxation time lower than 0.5 ( Chapman-Enskog for navier stokes suggestes us ) . We extract the viscosity from reynolds definition after adjusting characteristic velocity ( lid velocity) and characteristic length (ny- the number of lattice nodes in y direction) . By obtaning kinematic viscosity relaxation time could be obtained by satisfying the mentioned criteria .

2- We already set D2Q9 lattice velocities so we do not need to change streaming section however collision distribution function should be modified to account the advection roles . ( Make attention how we change this par)

ezafe kardane formule vase ghesmate collision

3- comparison were made by "Ghia.et al" . we used the u velocity along the horizontal line and V velocity along vertical line at the middle of the cavity .We got these data from " Gerris Flow Solver " site .

Further efforts :

1- Try to change the Reynolds number to see what would be the vortex behavior . 2- Try to change the number of nodes in Lattice . what is the relation between lattice nodes and the length of the cavity? 3- Try to change the lid driven position ( for example change the left side from no-slip to moving boundary) 4- Try to creat cavity with unequal sides' length . Lets the length of the cavity is twice of the width .


In [166]:
nx=100 # the number of nodes in x direction lattice direction 
ny=100 # the number of nodes in y direction lattice direction 
Re=1000 # rayleigh number
uo=0.1
rhoo=1.0
visco=uo*ny/Re
D=3.0 # the dimension of the problem 
omega=1./(0.5+(visco*D))     #Chapman-Enskog  dt=1 and dx=1  relaxation flow
mstep=100000 # the number of time step
k=9 # k=0,1,2,3,4,5,6,8,9

In [167]:
w=numpy.ones(k)      # witghting factor
cx=numpy.ones(k) 
cy=numpy.ones(k) 
rho=numpy.ones((ny+1,nx+1) )   # density matrix
u=numpy.ones((ny+1,nx+1) )  
v=numpy.ones((ny+1,nx+1) )  
f= numpy.ones((k, ny+1,nx+1))   # distribution function
feq= numpy.ones((k, ny+1,nx+1))
t1=numpy.ones((ny+1,nx+1) )   
t2=numpy.ones((ny+1,nx+1) )  
rhon=numpy.ones(nx+1)

In [168]:
y1= numpy.zeros(17)
u1= numpy.zeros(17)

y1[0]=-0.49933  ;u1[0]= -0.000882
y1[1]=-0.444335 ;u1[1]= -0.181701
y1[2]=-0.43629  ;u1[2]= -0.201989
y1[3]=-0.428914 ;u1[3]= -0.222276
y1[4]=-0.397406 ;u1[4]= -0.297251
y1[5]=-0.327052 ;u1[5]= -0.383699
y1[6]=-0.217948 ;u1[6]= -0.27788
y1[7]=-0.046595 ;u1[7]= -0.106804
y1[8]=0.001598  ;u1[8]= -0.060949
y1[9]=0.118733  ;u1[9]= 0.057217
y1[10]=0.235193 ;u1[10]= 0.186849
y1[11]=0.352315 ;u1[11]= 0.333239
y1[12]=0.45404  ;u1[12]= 0.466401
y1[13]=0.461386 ;u1[13]= 0.511382
y1[14]=0.469392 ;u1[14]= 0.574884
y1[15]=0.476719 ;u1[15]= 0.659554
y1[16]=0.5      ;u1[16]= 0.999118

In [169]:
x1= numpy.zeros(17)
v1= numpy.zeros(17)

x1[0]=-0.5 ;v1[0]=0.00069404
x1[1]=-0.43768  ;v1[1]=0.275621
x1[2]=-0.429602 ;v1[2]=0.290847
x1[3]=-0.421523 ;v1[3]=0.303994
x1[4]=-0.406521 ;v1[4]=0.326826
x1[5]=-0.343624 ;v1[5]=0.371038
x1[6]=-0.273803 ;v1[6]=0.330015
x1[7]=-0.265724 ;v1[7]=0.32307
x1[8]=-0.000289 ;v1[8]=0.0252893
x1[9]=0.304962  ;v1[9]=-0.318994
x1[10]=0.359781 ;v1[10]=-0.427191
x1[11]=0.40652  ;v1[11]=-0.515279
x1[12]=0.445182 ;v1[12]=-0.392034
x1[13]=0.45326  ;v1[13]=-0.336623
x1[14]=0.461339 ;v1[14]=-0.277749
x1[15]=0.46884  ;v1[15]=-0.214023
x1[16]=0.5      ;v1[16]=-6.20706e-17

In [170]:
strf=numpy.ones((ny+1,nx+1) )   # streamfunction 
##================ Initial boundary condition
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.

cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0

cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0

In [171]:
##================== Initial value

rho[0:ny+1,0:nx+1]=rhoo
v[0:ny+1,0:nx+1]=0.0
u[0:ny+1,0:nx+1]=0.0

u[ny,0:nx+1]=uo



for i in range(nx+1):
  for j in range(ny+1):
   for l in range (k): #k=0,1,2,3,4      
    f[l,j,i]=w[l]*rho[j,i]

In [172]:
##   Main loop  : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
 #print n ,u[nx/2,ny/2]
 if (n%10000==0):
  print n
# collision process fluid flow  
 
  
 t1[0:ny+1,0:nx+1]=u[0:ny+1,0:nx+1]*u[0:ny+1,0:nx+1]+v[0:ny+1,0:nx+1]*v[0:ny+1,0:nx+1]
 for l in range (k):
  t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
  feq[l,0:ny+1,0:nx+1]=w[l]*rho[0:ny+1,0:nx+1]*(  1.+ 3.*t2[0:ny+1,0:nx+1]+4.5*t2[0:ny+1,0:nx+1]**2-1.5*t1[0:ny+1,0:nx+1]   )    
  f[l,0:ny+1,0:nx+1]=(1.-omega)*f[l,0:ny+1,0:nx+1]+omega*feq[l,0:ny+1,0:nx+1] 
  
 f[1,0:ny+1,nx:0:-1]=f[1,0:ny+1,nx-1::-1]
 f[3,0:ny+1,0:nx]=f[3,0:ny+1,1:nx+1]
 f[2,ny:0:-1,0:nx+1]=f[2,ny-1::-1,0:nx+1]
 f[5,ny:0:-1,nx:0:-1]=f[5,ny-1::-1,nx-1::-1]
 f[6,ny:0:-1,0:nx]=f[6,ny-1::-1,1:nx+1] 
 f[4,0:ny,0:nx+1]=f[4,1:ny+1,0:nx+1]
 f[8,0:ny,nx:0:-1]=f[8,1:ny+1,nx-1::-1]
 f[7,0:ny,0:nx]=f[7,1:ny+1,1:nx+1]
 
 ##  =============================
                #left Boundary- Bounce back
 f[1,0:ny+1,0]=f[3,0:ny+1,0]
 f[5,0:ny+1,0]=f[7,0:ny+1,0]
 f[8,0:ny+1,0]=f[6,0:ny+1,0]
  
                 #right Boundary-bounce back
 f[3,0:ny+1,nx]=f[1,0:ny+1,nx]
 f[7,0:ny+1,nx]=f[5,0:ny+1,nx]
 f[6,0:ny+1,nx]=f[8,0:ny+1,nx]
  
               # top Boundary -bounce back
 rhon[0:nx+1]=f[0,ny,0:nx+1]+f[1,ny,0:nx+1]+f[3,ny,0:nx+1]+2*(f[2,ny,0:nx+1]+f[6,ny,0:nx+1]+f[5,ny,0:nx+1])               
 f[4,ny,0:nx+1]=f[2,ny,0:nx+1]
 f[7,ny,0:nx+1]=f[5,ny,0:nx+1]-(rhon[0:nx+1]*uo/6.0)
 f[8,ny,0:nx+1]=f[6,ny,0:nx+1]+(rhon[0:nx+1]*uo/6.0)
  

                 # bottom  Boundary -bunce back
 f[2,0,0:nx+1]=f[4,0,0:nx+1]
 f[5,0,0:nx+1]=f[7,0,0:nx+1]
 f[6,0,0:nx+1]=f[8,0,0:nx+1]
 
 for i in range(nx+1):
  for j in range(ny+1):
   sumr=0.0
   for l in range (k):
    sumr=sumr+f[l,j,i]
   rho[j,i]=sumr
   
 for i in range(1,nx):
  for j in range(1,ny):
   usum=0.0
   vsum=0.0
   for l in range (k):
    usum=usum+f[l,j,i]*cx[l]
    vsum=vsum+f[l,j,i]*cy[l]
   u[j,i]=usum/rho[j,i]
   v[j,i]=vsum/rho[j,i]
   
 
 
 u[:,0]=0.
 u[:,nx]=0.
 u[0,:]=0.
 u[ny,:]=uo
 
 v[:,0]=0.
 v[:,nx]=0.
 v[0,:]=0.
 v[ny,:]=0.


0
10000
20000
30000
40000
50000
60000
70000
80000
90000

In [173]:
strf[0,0]=0.0
for i in range(nx+1):
 rhoav=0.5*(rho[0,i-1]+rho[0,i])
 if(i>0) :
  strf[0,i]=strf[0,i-1]-rhoav*0.5*(v[0,i-1]+v[0,i] )
 for j in range(1,ny+1):
  rhom=0.5*(rho[j,i]+rho[j-1,i]) 
  strf[j,i]=strf[j-1,i]+rhom*0.5*(u[j,i]+u[j-1,i]) 
  
x=numpy.linspace(0,nx,nx+1)
y=numpy.linspace(0,ny,ny+1)

In [174]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,strf,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();



In [190]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contour(x,y,strf,900)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();



In [175]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,u,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();



In [176]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,v,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();



In [177]:
x1=(x1*ny)+(ny/2)
y1=(y1*ny)+(ny/2)

vv=v[ny/2,0:nx+1]/uo
uu=u[0:ny+1,nx/2]/uo

plt.plot(x1,v1, 'r*',label=' Ghia et.al'); 
plt.plot(x,vv, 'b-',label=' LBM '); 
plt.legend();



In [178]:
plt.plot(y1,u1, 'r*',label=' Ghia et.al'); 
plt.plot(y,uu, 'b-',label=' LBM '); 
plt.legend();